G D A R
GULF DATA, ASSESSMENT, AND REVIEW

GDAR 01
Gulf of Mexico Blue Crab

Model Code Supplement

Prepared by
S. VanderKooy (editor)

Gulf States Marine Fisheries Commission
2404 Government Street
Ocean Springs, Mississippi, 39564

Appendix A.1   ADMB Source Code for Two-stage Model (CMSA.tpl)

//###################################################################
//###################################################################

// Adapted from Chesapeake Bay 2010 blue crab assessment 
//By: M. Wilberg 2/28/2011
//
//Adapted by W Cooper 2012 
//Major changes:
//	1) Pulled out sex-specific due to limited biostatistical sampling of landings in Gulf 
//	2) Added stage-specific mortality, including in ref points
//	3) Added environmental influences on S-R and stage-specific M 
// 4) Added retrospective analyses and projections

//###################################################################
//###################################################################

TOP_OF_MAIN_SECTION
//increase number of estimated parameters
  gradient_structure::set_NUM_DEPENDENT_VARIABLES(2000);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(2000400);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(10000000);
  arrmblsize = 10000000;

//###################################################################
//###################################################################
DATA_SECTION
  init_adstring DataFile;
  !! ad_comm::change_datafile_name(DataFile);

  init_int testing	//toggle to turn on/off console output for testing
  init_int fyear   //first year of the model run
  init_int lyear   //last year of the model run
  init_int retroYears
  init_int projYears

  int timeSteps  //number of time steps in a year
  int stages	//number of stages/ages to model
  !! timeSteps=1;
  !! stages=2;

  int retroSteps
  int projSteps
  !! retroSteps=retroYears*timeSteps;
  !! projSteps=projYears*timeSteps;

  //Because coded to be multiple time steps per year, need to define an indexing scheme that isn't year-based
  int mTimeSteps	// total model time steps for population dynamics
  int startIndex	//starting index for the model
  int endIndex		//ending index for the model
  !!mTimeSteps = (lyear-fyear+1)*timeSteps; //calculate the time steps in the model
  !!startIndex=1000;		//set a value, can be anything > the lag of environment time series 

  //substract off the retrospective period
  !!endIndex=startIndex+mTimeSteps-1-retroSteps;  //substract 1 since startIndex is first step

  //#################################
  //#Catch Data
  //#################################
  init_int ftcyear  //first year of total catch
  init_int ltcyear  //last year of total catch
  int cStartIndex
  int cTimeSteps  //catch time steps based on catch years
  !!cTimeSteps = (ltcyear-ftcyear+1)*timeSteps;
  !!cStartIndex = (ftcyear-fyear)*timeSteps+startIndex;
  init_vector com_TC_obs(cStartIndex,cStartIndex+cTimeSteps-1)  //total catch
  init_number C_cv					//catch data standard deviation
  init_int effFlag		//flag to use effort time series
  init_vector com_Eff_obs(cStartIndex,cStartIndex+cTimeSteps-1)  //total catch

  //#################################
  //#Survey data
  //#################################
  //Adults
  init_int numAdSurv
  init_int fsayear  //first year of adult surveys
  init_int lsayear  //last year of adult surveys
  int adTimeSteps
  int adStartIndex
  !!adTimeSteps=(lsayear-fsayear+1)*timeSteps;
  !!adStartIndex = (fsayear-fyear)*timeSteps+startIndex;
  init_matrix ad_survey_obs(1,numAdSurv,adStartIndex,adStartIndex+adTimeSteps-1)  //Adult survey CPUE
  init_matrix ad_survey_cv(1,numAdSurv,adStartIndex,adStartIndex+adTimeSteps-1)  	//Survey CVs for adults
  init_vector sa_time(1,numAdSurv)		//survey time

  //Recruits
  init_int numRecSurv
  init_int fsryear
  init_int lsryear
  int recTimeSteps
  int recStartIndex
  !!recTimeSteps=(lsryear-fsryear+1)*timeSteps;
  !!recStartIndex = (fsryear-fyear)*timeSteps+startIndex;
  init_matrix re_survey_obs(1,numRecSurv,recStartIndex,recStartIndex+recTimeSteps-1)  //Recruit survey CPUE
  init_matrix re_survey_cv(1,numRecSurv,recStartIndex,recStartIndex+recTimeSteps-1)  	//survey SDs for recruits
  init_vector sr_time(1,numRecSurv)

  //#################################
  //#Fishery params
  //#################################
  init_number p_rec  //Proportion of recreational harvest per region
  init_number p_under  //Proportion of harvest underreporting per region
  init_number maxF	//Max F for F_pen calculation
  init_number maxM	//Max M for F_pen calculation

  //#################################
  //#Adult Z estimates as prior
  //#################################
  init_number aveZ
  init_number Z_cv

  //#################################
  //#Life History params
  //#################################
  init_number sratio					//Sex ratio
  init_vector M(1,stages)			//mortality at age for each stage (e.g., for CS: recruits, post-recruits)
  vector Myr(1,stages)			//M rate on per year basis (for ref points)
  init_vector pSpawn(1,stages)  	//proportion of females spawning in each season for differential spawning throughout the year
  init_number sp_time  //proportion of the time step before spawning occurs
  init_int SRSwitch				//switch for recruit function formulation; 1=bev holt, 2=ricker

 LOCAL_CALCS
//convert M to the appropriate time frame from per year basis
  Myr=M;
  M=M/timeSteps;
 END_CALCS

  //#################################
  //#Environmental time series params/data
  //#################################
  init_int numEnvTS
  init_int feyear   //first year of the model
  init_int leyear   //last year of the model
  int eTimeSteps
  int eStartIndex		//this should be less than the startIndex
  !!eTimeSteps = (leyear-feyear+1)*timeSteps;
  !!eStartIndex = (feyear-fyear)*timeSteps+startIndex;
  init_matrix envObs(1,numEnvTS,eStartIndex,eStartIndex+eTimeSteps-1)	//environmental time series (regions, season,timesteps)
  init_vector env_cv(1,numEnvTS)	//environmental time series (regions, season,timesteps)
  init_int envRecTS					//time series # that influences recruitment
  init_int envRecLag				//lag of recruitment influence
  init_vector envMTS(1,stages)					//time series # that influences mortality
  init_vector envMLag(1,stages)						//lag in mortality influence
  matrix env(1,numEnvTS,eStartIndex,endIndex+projSteps+1) //add on one for the forward stepping recruitment

  //#################################
  //#Projections
  //#################################
  init_matrix envProj(1,numEnvTS,endIndex+1,endIndex+projSteps+1)
  init_vector effProj(endIndex+1,endIndex+projSteps+1) //START with terminal year since non-estimable F

  //#################################
  //#Parameters (initial val, min, max, phase)
  //#################################
  
  init_number init_NIni
  init_vector init_NParams(1,3)
  init_number init_RIni
  init_vector init_RParams(1,3)

  //F params
  init_number F_qIni
  init_vector F_qParams(1,3)
  init_number F_devIni
  init_vector F_devParams(1,3)
  init_number eff_cvIni
  init_vector eff_cvParams(1,3)

  //recruitment  params
  init_number rec_devIni
  init_vector rec_devParams(1,3)
  init_number rec_cvIni
  init_vector rec_cvParams(1,3)

  init_number S0Ini
  init_vector S0Params(1,3)
  init_number steepIni
  init_vector steepParams(1,3)

  init_number sr_beta_envIni
  init_vector sr_beta_envParams(1,3)
  init_vector M_beta_envIni(1,stages)
  init_vector M_beta_envParams(1,3)
  init_number M_cvIni
  init_vector M_cvParams(1,3)
  init_vector selIni(1,stages)
  init_vector selParams(1,3)

  //#################################
  //#Likelihood Weights
  //#################################
  
  init_number com_lambda		//survey weight
  init_vector sa_lambda(1,numAdSurv)		//survey weight
  init_vector sr_lambda(1,numRecSurv)
  init_number recDev_lambda		//survey weight
  init_number effort_lambda		//survey weight
  init_number aveZ_lambda		//survey weight

  //#################################
  //#Additional param control flags not addressed in data section
  //#################################
  init_number biasAdj					//Adjustment multiplier for bias correction factor

  //#################################
  //#Reference point calcs
  //#################################
  //number for reference point explorations
  init_number Fval_init  							//lowest value of F used in SPR calcs
  init_number Fval_max   						//highest value of F used in SPR calcs
  init_number Fval_inc   							//increment for F
  int Fval_num
  !!Fval_num=(Fval_max-Fval_init)/Fval_inc+1;
  init_int nspr          //number of values F-%SPR will be calculated at
  init_vector SPR_targ(1,nspr)  //Values of SPR for Fval reference point calculations

  //#################################
  //#EOF Test
  //#################################
  init_int test			//check that data read in appropriately

  //#################################
  //#Additional Variables
  //#################################
  
  //Total harvest including recreational
  vector TC_obs(cStartIndex,cStartIndex+cTimeSteps-1)  //total catch

  //Variances for data sets
  number C_var  //variance of catch
  number Z_var  //variance of catch
  matrix ad_survey_var(1,numAdSurv,adStartIndex,adStartIndex+adTimeSteps-1)  //variances for adult surveys
  matrix re_survey_var(1,numRecSurv,recStartIndex,recStartIndex+recTimeSteps-1)  //variances for recruitment surveys

  //Define index varaibles
  int y  //index variable for time step
  int s  //index variable for season
  int r  //index variable for region
  int i  //index variable 
  number year  //for report section
  int ispr
  int iter
  int iterMCMC
  !!iterMCMC=0;
  int index

 LOCAL_CALCS

  if (SRSwitch==2) steepParams(2)=5.0;	//make sure to bound steepness appropriately for Ricker

  C_var=log(C_cv*C_cv+1);  //variance of catch
  Z_var=log(Z_cv*Z_cv+1);  //variance of Z

  //commercial + recreational catch+prop. underreported 
  TC_obs=com_TC_obs*(1.+p_rec+p_under);

  //Calculate variances from CVs
  for (i=1; i<=numAdSurv; i++){
    for (y=adStartIndex; y<=adStartIndex+adTimeSteps-1; y++){
      ad_survey_var(i,y)=log(ad_survey_cv(i,y)*ad_survey_cv(i,y)+1);  //variances for adult surveys
    }
  }
  for (i=1; i<=numRecSurv; i++) {
    for (y=recStartIndex; y<=recStartIndex+recTimeSteps-1; y++){
      re_survey_var(i,y)=log(re_survey_cv(i,y)*re_survey_cv(i,y)+1);  //variances for recruitment surveys
    }
  }

  for (i=1; i<=numEnvTS; i++){
    for (y=eStartIndex; y<=endIndex; y++){
      env(i,y)=envObs(i,y);
    }
    for (y=endIndex+1; y<=endIndex+projSteps; y++){
      env(i,y)=envProj(i,y);
    }
  }

  if(test!=12345)  //check to make sure end of file number is correct
  {
    //if not correct, output the data and exit.
    cout << "Data not reading properly" << endl;
    cout << "Commercial\t" << com_TC_obs << endl;
    cout << "adults\t" <<ad_survey_obs << endl;
    cout << "recruits\t" <<re_survey_obs << endl;
    cout << "environment\t" <<env << endl;
    cout << "max F\t" <<maxF << endl;
    cout << "max M\t" <<maxM << endl;
    cout << "envMTS\t" <<envMTS << endl;
    cout << "env Lag\t" <<envMLag << endl;
    cout << "ini N\t" <<init_NIni << endl;
    cout << "steepIni\t" <<steepIni << endl;
    cout << "sel params\t" <<selParams<< endl;
    cout << "EOF test: " << test << endl;
    exit(1);
  }

 END_CALCS

//###################################################################
//###################################################################
PARAMETER_SECTION

//Copy the Parameter estimates to double vals so can put as bounds,phases below
 LOCAL_CALCS
//template: double xxxMin=log(xxxParams(1)); double xxxMax=log(xxxParams(2)); double xxxPhase=log(xxxParams(2));
  double log_init_NMin=log(init_NParams(1)); 			double log_init_NMax=log(init_NParams(2)); 			double init_NPhase=init_NParams(3);
  double log_init_RMin=log(init_RParams(1)); 			double log_init_RMax=log(init_RParams(2)); 			double init_RPhase=init_RParams(3);
  double log_F_qMin=log(F_qParams(1)); 			double log_F_qMax=log(F_qParams(2)); 		double F_qPhase=F_qParams(3);
  double log_F_devMin=log(F_devParams(1)); 			double log_F_devMax=log(F_devParams(2)); 			double F_devPhase=F_devParams(3);
  double log_eff_cvMin=log(eff_cvParams(1)); 		double log_eff_cvMax=log(eff_cvParams(2)); 			double eff_cvPhase=eff_cvParams(3);
  double log_rec_devMin=log(rec_devParams(1)); 	double log_rec_devMax=log(rec_devParams(2));		double rec_devPhase=rec_devParams(3);
  double log_rec_cvMin=log(rec_cvParams(1)); 		double log_rec_cvMax=log(rec_cvParams(2)); 			double rec_cvPhase=rec_cvParams(3);
  double log_S0Min=log(S0Params(1)); 					double log_S0Max=log(S0Params(2)); 						double S0Phase=S0Params(3);
  double log_steepMin=log(steepParams(1)); 			double log_steepMax=log(steepParams(2)); 				double steepPhase=steepParams(3);
  double sr_beta_envMin=sr_beta_envParams(1); 		double sr_beta_envMax=sr_beta_envParams(2); 		double sr_beta_envPhase=sr_beta_envParams(3);
  double M_beta_envMin=M_beta_envParams(1); 		double M_beta_envMax=M_beta_envParams(2); 		double M_beta_envPhase=M_beta_envParams(3);
  double log_M_cvMin=log(M_cvParams(1)); 		double log_M_cvMax=log(M_cvParams(2)); 			double M_cvPhase=M_cvParams(3);
  double log_selMin=log(selParams(1)); 					double log_selMax=log(selParams(2)); 					double selPhase=selParams(3);

 END_CALCS

  //initial R and N
  //template: (Min,Max,Phase)
  init_bounded_number log_init_N(log_init_NMin,log_init_NMax,init_NPhase)
  init_bounded_number log_init_R(log_init_RMin,log_init_RMax,init_RPhase)

  //Fishing mortality for each year
  init_bounded_number log_F_q(log_F_qMin,log_F_qMax,F_qPhase)
  init_bounded_vector log_F_dev(startIndex+1,endIndex-1,log_F_devMin,log_F_devMax,F_devPhase) //don't bother with terminal year deviation
  init_bounded_number log_eff_cv(log_eff_cvMin,log_eff_cvMax,eff_cvPhase)

  //Recruitment params
  init_bounded_dev_vector log_rec_dev(startIndex+1,endIndex,log_rec_devMin,log_rec_devMax,rec_devPhase)
  init_bounded_number log_rec_cv(log_rec_cvMin,log_rec_cvMax,rec_cvPhase)

  //Stock-recruitment parameters
  init_bounded_number log_S0(log_S0Min,log_S0Max,S0Phase)
  init_bounded_number log_steep(log_steepMin,log_steepMax,steepPhase)

  //S-R environmental link parameter
  init_bounded_number sr_beta_env(sr_beta_envMin,sr_beta_envMax,sr_beta_envPhase)

  //Adult environmental link parameter
  init_bounded_vector M_beta_env(1,stages,M_beta_envMin,M_beta_envMax,M_beta_envPhase)
  init_bounded_number log_M_cv(log_M_cvMin,log_M_cvMax,M_cvPhase)

  //Vulnerability at each stage
  init_bounded_vector log_sel(1,stages,log_selMin,log_selMax,selPhase)

  //############### Derived parameters ###############//
  
  //sdreport_matrix  for some of these
  sdreport_vector N(startIndex,endIndex+1+projSteps)   //abundance
  sdreport_vector R(startIndex,endIndex+1+projSteps)   //recruitment
  vector SP(startIndex,endIndex+projSteps) //number of spawners
  vector TC(startIndex,endIndex+projSteps)      //total catch
  vector u(startIndex,endIndex+projSteps)      //total exploitation rate
  sdreport_vector F(startIndex,endIndex+projSteps)   //fishing mortality rate
  vector effort(startIndex,endIndex+projSteps)   //fishing mortality rate
  vector sel(1,stages)				//selectivity (partial recruitment) of recruits to the fishery
  matrix Mt(1,stages,startIndex,endIndex+projSteps)		//M at time t; don't do for years to account for seasonal M
  vector Z(startIndex,endIndex+projSteps)		//total Z

  matrix ad_survey_est(1,numAdSurv,startIndex,endIndex)  //estimated adult survey indices
  matrix re_survey_est(1,numRecSurv,startIndex,endIndex)  //estimated recruitment survey indices
  vector qa(1,numAdSurv)   //catchability for adult surveys 
  vector qr(1,numRecSurv)   //catchability for recruitment surveys

  number S0
  number steep
  number R0
  number A0
  number alpha			//Alpha of the S-R relationship
  number beta			//Beta of the S-R relationship
  number rec_var		//variance for recruitment deviations
  number eff_var			//variance for effort residuals 
  number M_var			//variance for F deviations

  //Derived variables for reference point calculations 
  number Fval							//F for SPR calculations
  vector SPR(1,Fval_num)				//spawners per recruit (NOT spawning potential ratio)
  number SPR0							//virgin SPR
  number SPR1							//spawner per recruit per year
  vector tSPR(startIndex,endIndex+projSteps) //transitional SPR
  vector NPR(1,Fval_num)			//numbers per reruit
  vector YPR(1,Fval_num)				//yield per recruit
  vector SPRatio(1,Fval_num)		//spawning potential ratio
  vector N_eq(1,Fval_num)			//equilibruim numbers
  vector R_eq(1,Fval_num)			//equilibrium recruitment
  vector C_eq(1,Fval_num)			//equilibrium catch
  vector C_eqSort(1,Fval_num)			//sorted equilibrium catch

  //Reference points  
  vector u0_eq(1,Fval_num)			//equilibrium explotation rate for age0+
  vector u1_eq(1,Fval_num)			//equilibrium explotation rate for age1+
  vector uAll_eq(1,Fval_num)			//equilibrium explotation rate for age0+
  vector FSPR_ref(1,nspr)				//F%SPR reference points
  vector SPRDiff(1,nspr)				//temporary array to check if at proper F for SPR calcs
  vector Fvec(1,Fval_num)			//equilibrium explotation rate for age0+

  likeprof_number MSY				//MSY estimate
  number OFL				//Overfishing Limit (Ncurrent*uMSY)
  number u0MSY			//exploitation rate at MSY for age 0
  number u1MSY			//exploitation rate at MSY for age 1+
  number uMSY			//exploitation rate at MSY for age 0+
  number FMSY			//F rate at MSY
  number RMSY		//equilibrium recruitment at msy
  number NMSY			//Number at MSY
  number FLim			//F Limit (target)
  number NLim			//N Limit (target)
  number cLim			//c used in calculation of targets
  number FCurr			//Current F (geometric mean of last 3 years of model run, not including terminal year)
  vector MCurr(1,stages)			//Current M (geometric mean of last 3 years of model run)
  number NCurr			//Current N (geometric mean of last 3 years of model run)
  number SPRCurr			//Current SPR using FCurr and MCurr
  number FMSYRatio	//Fcurr/FMSY
  number NMSYRatio	//Ncurr/NMSY
  number UMSYRatio	//Ncurr/NMSY
  number FFLimRatio	//Fcurr/FLimit
  number NNLimRatio	//Ncurr/NLim
  number termToMSY  //for F calcs 

  //variables for likelihood function
  vector Lsr(1,numRecSurv)  //likelihood components for recruit surveys
  vector Lsa(1,numAdSurv)  //likelihood components for adult surveys
  number Lc        //likelihood components for catch time series 
  number Lz        //likelihood components for adult Z estimate (anchors F/N0/R0)  
  number Lrdev           //likelihood for recruitment deviations
  number Leff           //likelihood for effort residuals

  number F_pen           //penalty for F above the max F
  number M_pen           //penalty for M above the max M
  objective_function_value negLL 	//negative log-liklihood

  //############### Starting parameter values ###############//
 LOCAL_CALCS

  log_init_N=log(init_NIni);
  log_init_R=log(init_RIni);
  log_F_q=log(F_qIni);
  log_F_dev=log(F_devIni);
  log_M_cv=log(M_cvIni);
  log_rec_cv=log(rec_cvIni);
  log_rec_dev=log(rec_devIni);
  log_S0=log(S0Ini);
  log_steep=log(steepIni);
  sr_beta_env=sr_beta_envIni;
  M_beta_env=M_beta_envIni;
  log_sel=log(selIni);

 END_CALCS

//###################################################################
//###################################################################
PROCEDURE_SECTION
  set_initial_conditions();
  if (testing==1)  cout << "End set_initial_conditions()" << endl;
  calculate_abundance_and_catch();
  if (testing==1)  cout << "End calculate_abundance_and_catch()" << endl;
  calculate_predicted_indices();
  if (testing==1)  cout << "End calculate_predicted_indices()" << endl;
  calculate_objective_function();
  if (testing==1)  cout << "End calculate_objective_function()" << endl;
  mcmc();
  if (testing==1) cout << "End mcmc()" << endl;

  if (testing==1) {
    calculate_tSPR();
    obs_pred();
    MSY_estimates();
    HPD_estimates();
    general_report();

    cout << "Procedure section completed first cycle, now exiting"<< endl;
    exit(1); //exit if in testing phase -- runs model at initial parameter values
  }

//############### Main Functions ###############

FUNCTION set_initial_conditions

  //convert parameters from the log scale
  S0=exp(log_S0);
  steep=exp(log_steep);
  negLL=0.0;
  M_var=log(exp(log_M_cv)*exp(log_M_cv)+1);
  rec_var=log(exp(log_rec_cv)*exp(log_rec_cv)+1);
  eff_var=log(exp(log_eff_cv)*exp(log_eff_cv)+1);
  sel=exp(log_sel);
  F_pen=0;
  M_pen=0;

  //######## S-R Params ########  
  
  //Calculate virgin SPR, including proportion of recruits spawning  
  A0=sratio*pSpawn(2)*exp(-(M(1)+sp_time*M(2)))/(1.-exp(-(M(2)))) + sratio*pSpawn(1)*exp(-(sp_time*M(1)));
  R0=S0/A0;

  if (SRSwitch==1)	{ //Beverton-Holt
    alpha = S0*(1-steep)/(4*steep*R0);
    beta = (5*steep-1)/(4*steep*R0);
  }

  if (SRSwitch==2)	{ //Ricker
    beta = log(5*steep)/(0.8*S0);
    alpha =(exp((5.*log(5.*steep))/4.)*R0)/S0;
  }

  //calculate reference points after setting S-R params so can get FMSY for F projections
  calculate_reference_points();

  //######## M ########  
  
  //compute the yearly M accounting for seasonal differences and environmental differences
  //leave this here to deal with seasonality
  for(y=startIndex; y<=endIndex+projSteps; y++) {

    //only apply deviation + bias correction if active
    Mt(1,y)=M(1);
    if (active(M_beta_env)) Mt(1,y)=M(1)*exp(M_beta_env(1)*env(envMTS(1),y-envMLag(1)))*exp(-0.5*M_var);

    if (y<=endIndex) {
      posfun(maxM-(timeSteps*Mt(1,y)),.000001,M_pen);
      negLL+=100.*M_pen;
    }

    Mt(2,y)=M(2);
    if (active(M_beta_env)) Mt(2,y)=M(2)*exp(M_beta_env(2)*env(envMTS(2),y-envMLag(2)))*exp(-0.5*M_var);

    if (y<=endIndex) {
      posfun(maxM-(timeSteps*Mt(2,y)),.000001,M_pen);
      negLL+=100.*M_pen;
    }
  }

  //######## F ########  
  
  if (effFlag==0) effort=1.0;
  else {
    //set up effort to average and replace all missing data with average
    double avg_effort=mean(com_Eff_obs);

    //for any year prior to effort data, set to avg of all other years
    for (i=startIndex; i<=endIndex; i++) effort(i)=com_Eff_obs(i);
    for (i=endIndex+1; i<=endIndex+projSteps; i++) effort(i)=avg_effort; //effort(endIndex)+effProj(i)*termToMSY; //deviation off the last year
    effort/=avg_effort; //scale to observed years and not including projected years
  }

  //If want to estimate ave. q instead of 1st year q: change F_dev to bounded_dev_vector and adjust here
  F(startIndex)=exp(log_F_q+log(effort(startIndex)));
  posfun(maxF-(timeSteps*F(startIndex)),.000001,F_pen);
  negLL+=100.*F_pen;

  for(y=startIndex+1; y<endIndex; y++) { //don't include terminal year estimate
    //Computed as F=q*Eff*exp(dev)
    F(y)=exp(log_F_q+log(effort(y))+log_F_dev(y));
    posfun(maxF-(timeSteps*F(y)),.000001,F_pen);
    negLL+=100.*F_pen;
  }

  //for terminal year, use estimated deviation from previous year to keep scaled together
  F(endIndex)=exp(log_F_q+log(effort(endIndex))+log_F_dev(endIndex-1));

  termToMSY=FMSY-F(endIndex); //effort range from terminal year to MSY; negative if FMSY<termF 

  //no F deviations on terminal year or projected years
  for(y=endIndex+1; y<=endIndex+projSteps; y++) {
    F(y)=(F(endIndex)+effProj(y)*termToMSY);
    effort(y)=exp(log(F(y)))/exp(log_F_q);
  }

  F=F/timeSteps;

  //######## Adult Z ########  
  
  for(y=startIndex; y<=endIndex+projSteps; y++) Z(y)=F(y)+Mt(2,y);

FUNCTION calculate_abundance_and_catch

  N(startIndex)=exp(log_init_N);
  R(startIndex)=exp(log_init_R);

  for(y=startIndex; y<=endIndex+projSteps; y++) {

    //spawners also include some animals that were recruits in the beginning of the year
    SP(y)=sratio*(N(y)*exp(-sp_time*(Mt(2,y)+sel(2)*F(y))))*pSpawn(2)
      + sratio*(R(y)*exp(-sp_time*(Mt(1,y)+sel(1)*F(y))))*pSpawn(1);

    if (SRSwitch==1)	{ //Beverton-Holt
      // don't use recruit deviations for projection years
      if (y<endIndex) R(y+1)=(SP(y)/(SP(y)*beta+alpha))*exp(sr_beta_env*env(envRecTS,y+1-envRecLag))*exp(log_rec_dev(y+1)-biasAdj*0.5*rec_var);
      else R(y+1)=(SP(y)/(SP(y)*beta+alpha))*exp(sr_beta_env*env(envRecTS,y+1-envRecLag));
    }

    if (SRSwitch==2)	{ //Ricker
      // don't use recruit deviations for projection years
      if (y<endIndex) R(y+1)=(alpha*SP(y)*exp(-beta*SP(y)))*exp(sr_beta_env*env(envRecTS,y+1-envRecLag))*exp(log_rec_dev(y+1)-biasAdj*0.5*rec_var);
      else R(y+1)=(alpha*SP(y)*exp(-beta*SP(y)))*exp(sr_beta_env*env(envRecTS,y+1-envRecLag));
    }

    //abundance for the next year
    N(y+1)=R(y)*exp(-(Mt(1,y)+sel(1)*F(y)))+N(y)*exp(-(Mt(2,y)+sel(2)*F(y)));

    //Baranov catch equation
    TC(y)=N(y)*((sel(2)*F(y))/(sel(2)*F(y)+Mt(2,y))*(1.-exp(-(sel(2)*F(y)+Mt(2,y)))))
      + R(y)*((sel(1)*F(y))/(sel(1)*F(y)+Mt(1,y))*(1.-exp(-(sel(1)*F(y)+Mt(1,y)))));

    u(y)=TC(y)/(R(y)*((1-exp(-sel(1)*F(y)))/(1-exp(-F(y))))+N(y));

  }

  //Calculate year-dependent F/N Reference Point components (i.e., ratios)
  NCurr=mfexp((log(N(endIndex))+log(N(endIndex-1))+log(N(endIndex-2)))/3);
  FCurr=mfexp((log(F(endIndex-1))+log(F(endIndex-2)))/2);
  FMSYRatio=FCurr/FMSY;
  NMSYRatio=NCurr/NMSY;
  UMSYRatio=mfexp((log(u(endIndex))+log(u(endIndex-1))+log(u(endIndex-2)))/3)/uMSY;
  OFL=uMSY*mfexp((log(N(endIndex))+log(N(endIndex-1))+log(N(endIndex-2)))/3);  //N is by region, so need to take mean if more

  cLim=max(1-M(2),0.5);
  NLim=cLim*NMSY;
  FLim=FMSY;
  if (NCurr <= NLim) FLim=(FMSY*NCurr)/(cLim*NMSY);
  FFLimRatio=FCurr/FLim;
  NNLimRatio=NCurr/NLim;

  MCurr(1)=mfexp((log(Mt(1,endIndex))+log(Mt(1,endIndex-1))+log(Mt(1,endIndex-2)))/3);
  MCurr(2)=mfexp((log(Mt(2,endIndex))+log(Mt(2,endIndex-1))+log(Mt(2,endIndex-2)))/3);

  SPR0=sratio*pSpawn(2)*exp(-(MCurr(1)+sp_time*MCurr(2)))/(1.-exp(-(MCurr(2))))
    + sratio*pSpawn(1)*exp(-(sp_time*MCurr(1)));
  SPR1=sratio*pSpawn(2)*exp(-(MCurr(1)+sel(1)*FCurr+sp_time*(MCurr(2)+sel(2)*FCurr)))/(1.-exp(-(MCurr(2)+sel(2)*FCurr)))
    + sratio*pSpawn(1)*exp(-sp_time*(MCurr(1)+sel(1)*FCurr));
  SPRCurr=SPR1/SPR0;

FUNCTION calculate_predicted_indices

  //########## Recruits #########
  for (i=1; i<=numRecSurv; i++){

    qr(i)=0.0;
    double counter=0.0;

    for(y=startIndex; y<=endIndex; y++) {
      if (y<recStartIndex) continue;

      if(re_survey_obs(i,y)!=-999.) { //check to make sure year is not missing
        if(!last_phase()) {
          //small constant added to recruitment in earlier stages to 
          //increase numerical stability
          //NOTE: this formulation assumes survey occurs before vulnerable to harvest
          qr(i)+=log(re_survey_obs(i,y))-log(R(y)*exp(-(sr_time(i)*Mt(1,y)))+.000001);
        }
        else { //small constant not included in last estimation stage
          qr(i)+=log(re_survey_obs(i,y))-log(R(y)*exp(-(sr_time(i)*Mt(1,y))));
        }
        counter++;
      }
    }
    //calculate geometric mean
    qr(i)=exp(qr(i)/counter);
    //Calculate predicted index of abundance
    //NOTE: this formulation assumes survey occurs before vulnerable to harvest
    for(y=startIndex; y<=endIndex; y++) {
      re_survey_est(i,y)=qr(i)*(R(y)*exp(-(sr_time(i)*Mt(1,y))));
    }
  }

  //########## Adults #########
  for (i=1; i<=numAdSurv; i++){
    //calculate catchability for each sex-index combination
    double counter=0.0;
    qa(i)=0.0;

    for(y=startIndex; y<=endIndex; y++) {
      if (y<adStartIndex) continue;

      if(ad_survey_obs(i,y)!=-999.) { //check to make sure year is not missing

        qa(i)+=log(ad_survey_obs(i,y))-log(N(y)*exp(-sa_time(i)*(Mt(2,y)+sel(2)*F(y))));
        counter++;
      }
    }
    //calculate geometric mean
    qa(i)=exp(qa(i)/counter);

    //Calculate each predicted index of abundance
    for(y=startIndex; y<=endIndex; y++) {
      ad_survey_est(i,y)=qa(i)*(N(y)*exp(-sa_time(i)*(Mt(2,y)+sel(2)*F(y))));
    }
  }

FUNCTION calculate_objective_function

  double pi=3.141593;

  //calculate adult survey likelihood component
  for (i=1; i<=numAdSurv; i++){
    Lsa(i)=0.0;
    for(y=startIndex; y<=endIndex; y++) {
      if (y<adStartIndex) continue;
      if(ad_survey_obs(i,y)!=-999.) { //check to make sure year is not missing -- some holes
        Lsa(i)+=0.5*log(2.*pi)+0.5*log(ad_survey_var(i,y))+log(ad_survey_obs(i,y))+square(log(ad_survey_obs(i,y)+.000001)-log(ad_survey_est(i,y)+.000001))/(2*ad_survey_var(i,y));
      }
    }
    Lsa(i)=sa_lambda(i)*Lsa(i);
  }

  //calculate recruit survey likelihood component
  for (i=1; i<=numRecSurv; i++){
    Lsr(i)=0.0;
    for(y=startIndex; y<=endIndex; y++) {
      if (y<recStartIndex) continue;
      if(re_survey_obs(i,y)!=-999.) { //check to make sure year is not missing
        Lsr(i)+=0.5*log(2.*pi)+0.5*log(re_survey_var(i,y))+log(re_survey_obs(i,y))+square(log(re_survey_obs(i,y)+.000001)-log(re_survey_est(i,y)+.000001))/(2*re_survey_var(i,y));
      }
    }
    Lsr(i)=sr_lambda(i)*Lsr(i);
  }

  //calculate total catch likelihood component
  Lc=0.0;
  for(y=startIndex; y<=endIndex; y++) {
    if (y<cStartIndex) continue;
    if(TC_obs(y)!=-999.) { //check to make sure year is not missing
      Lc+=square(log(TC_obs(y)+.000001)-log(TC(y)+.000001));
    }
  }
  Lc=com_lambda*(0.5*log(2.*pi)*size_count(TC_obs(startIndex,endIndex))+0.5*log(C_var)*size_count(TC_obs(startIndex,endIndex))+sum(log(TC_obs(startIndex,endIndex)))+0.5*Lc/C_var);

  //calculate likelihood component for recruitment deviations
  Lrdev=recDev_lambda*(0.5*log(2.*pi)*size_count(log_rec_dev)+0.5*log(rec_var)*size_count(log_rec_dev)+sum(log_rec_dev)+0.5*norm2(log_rec_dev)/rec_var);

  //calculate likelihood component for effort residuals if effort time series is included
  Leff=0.0;
  if (effFlag==1) {
    for(y=startIndex; y<endIndex; y++) {
      Leff+=0.5*log(2.*pi)+0.5*log(eff_var)+log(effort(y))+0.5*square(log(effort(y))-(log(F(y))-log_F_q))/eff_var;
    }
    Leff=effort_lambda*Leff;
  }

  //calculate likelihood component for total Z of adults as prior, read from independent Z estimate
  Lz=aveZ_lambda*(0.5*log(2.*pi)+0.5*log(Z_var)+log(aveZ)+0.5*square(log(aveZ)-log(mean(Z(startIndex,endIndex-1))))/Z_var);

  negLL+=sum(Lsa)+sum(Lsr)+Lc+Lrdev+Leff+Lz;

FUNCTION calculate_reference_points

  //Reference point variables
  MSY=0.0;
  u1MSY=0.0;
  u0MSY=0.0;
  uMSY=0.0;
  i=0;
  OFL=0.0;
  SPRDiff=1e10;

  //With recruit spawners 
  SPR0=sratio*pSpawn(2)*(exp(-(Myr(1)+sp_time*Myr(2)))/(1.-exp(-(Myr(2)))))
    + sratio*pSpawn(1)*(exp(-((sp_time*Myr(1)))));

  Fval=Fval_init;
  for(i=1; i<=Fval_num; i++)
  {
    Fvec(i)=Fval; //record the F values

    SPR(i)=sratio*pSpawn(2)*(exp(-(Myr(1)+sel(1)*Fval+sp_time*(Myr(2)+sel(2)*Fval)))/(1.-exp(-(Myr(2)+sel(2)*Fval))))
      + sratio*pSpawn(1)*(exp(-(sp_time*(Myr(1)+sel(1)*Fval))));

    NPR(i)=exp(-(Myr(1)+sel(1)*Fval))/(1.-exp(-(Myr(2)+sel(2)*Fval)));

    YPR(i)=(sel(1)*Fval)/(sel(1)*Fval+Myr(1))*(1.-exp(-(sel(1)*Fval+Myr(1))))
      + ((sel(2)*Fval)/(sel(2)*Fval+Myr(2))*(1.-exp(-(sel(2)*Fval+Myr(2)))))*NPR(i);

    if (SRSwitch==1)	 R_eq(i)=(SPR(i)-alpha)/(SPR(i)*beta);
    if (SRSwitch==2)	R_eq(i)=(log(alpha)+log(SPR(i)))/(beta*SPR(i));

    N_eq(i) = NPR(i)*R_eq(i);
    C_eq(i)=YPR(i)*R_eq(i);

    //calculate exploitation rate
    //age 0+
    u0_eq(i)=(sel(1)*Fval)/(sel(1)*Fval+Myr(1))*(1.-exp(-(sel(1)*Fval+Myr(1))));
    //age 1+
    u1_eq(i)=(sel(2)*Fval)/(sel(2)*Fval+Myr(2))*(1.-exp(-(sel(2)*Fval+Myr(2))));
    //all ages
    if (i>1) uAll_eq(i)=C_eq(i)/(N_eq(i)+R_eq(i)*((1-exp(-sel(1)*Fval))/(1-exp(-Fval))));

    //MSY
    if (C_eq(i)>MSY) {
      MSY=C_eq(i);
      FMSY=Fval;
      NMSY=N_eq(i);
      RMSY=R_eq(i);
      u0MSY=u0_eq(i);
      u1MSY=u1_eq(i);
      uMSY=uAll_eq(i);
    }

    //loop through SPR targets and see if at the correct F for each target
    for (ispr=1; ispr<=nspr; ispr++){
      if (square(SPR(i)/SPR0-SPR_targ(ispr)) < SPRDiff(ispr)) {
        SPRDiff(ispr)=square(SPR(i)/SPR0-SPR_targ(ispr));
        FSPR_ref(ispr)=Fval;
      }
    }

    //increment the female F for the SPR
    Fval+=Fval_inc;
  }

//############### Reporting functions ###############//

FUNCTION mcmc
  //Code to write results of MCMC to file so we can access the chains
  if(mceval_phase()) {
    //Define output file stream for MCMC results
    if(iterMCMC==0) {
      ofstream mcmcout("cmsa_refs.mcmc");
      mcmcout <<"MSY\t"<<"FMSY\t"<< "NMSY\t"<<"uMSY\t"<<"FLim\t"<<"NLim\t"<<"FMSYRatio\t"<<"NMSYRatio\t"<<"UMSYRatio\t"<<"FFLimRatio\t"<<"NNLimRatio\t"<<"SPRCurrent"<<endl;

      //print out yearly F and N
      ofstream mcmcout2("cmsa_yearly.mcmc");

      year=fyear;
      for(y=startIndex;y<=endIndex;y++){
        mcmcout2 <<"N"<<year<<"\t";
        year=year+1.0/timeSteps;
      }

      year=fyear;
      for(y=startIndex;y<=endIndex;y++){
        mcmcout2 <<"R"<<year<<"\t";
        year=year+1.0/timeSteps;
      }

      year=fyear;
      for(y=startIndex;y<=endIndex;y++){
        if( y<endIndex) mcmcout2 <<"F"<<year<<"\t";
        else mcmcout2 <<"F"<<year << endl;
        year=year+1.0/timeSteps;
      }

      ofstream mcmcout3("cmsa_pars.mcmc");
      mcmcout3 <<"N0\t"<<"R0\t"<< "Fq\t"<<"S0\t"<<"h"<<endl;

      iterMCMC++;
    }

    ofstream mcmcout("cmsa_refs.mcmc",ios::app);
    mcmcout <<MSY<<"\t"<<FMSY<<"\t"<<NMSY<<"\t"<<uMSY<<"\t"<<FLim<<"\t"<<NLim<<"\t"<<FMSYRatio<<"\t"<<NMSYRatio<<"\t"<<UMSYRatio<<"\t"<<FFLimRatio<<"\t"<<NNLimRatio<<"\t"<<SPRCurr<<endl;

    //print out yearly F and N
    ofstream mcmcout2("cmsa_yearly.mcmc",ios::app);

    for(y=startIndex;y<=endIndex;y++){
      mcmcout2 <<N(y) << "\t";
    }

    for(y=startIndex;y<=endIndex;y++){
      mcmcout2 <<R(y) << "\t";
    }

    for(y=startIndex;y<=endIndex;y++){
      if (y<endIndex) mcmcout2 <<F(y) << "\t";
      else mcmcout2 <<F(y) <<endl;
    }

    //print out yearly F and N
    ofstream mcmcout3("cmsa_pars.mcmc",ios::app);
    mcmcout3 <<exp(log_init_N)<<"\t"<<exp(log_init_R) <<"\t"<<exp(log_F_q)<<"\t"<<exp(log_S0)<<"\t"<<exp(log_steep) << endl;

  }

FUNCTION calculate_tSPR

  ofstream ofs_tSPR("tSPR.dat");
  ofs_tSPR<< "year\ttSPR" << endl;
  year=fyear;

  for(y=startIndex;y<=endIndex+projSteps;y++){

    SPR0=sratio*pSpawn(2)*exp(-(Mt(1,y)+sp_time*Mt(2,y)))/(1.-exp(-(Mt(2,y))))
      + sratio*pSpawn(1)*exp(-(sp_time*Mt(1,y)));
    SPR1=sratio*pSpawn(2)*exp(-(Mt(1,y)+sel(1)*F(y)+sp_time*(Mt(2,y)+sel(2)*F(y))))/(1.-exp(-(Mt(2,y)+sel(2)*F(y))))
      + sratio*pSpawn(1)*exp(-sp_time*(Mt(1,y)+sel(1)*F(y)));

    tSPR(y)=SPR1/SPR0;
    ofs_tSPR << year << "\t" << tSPR(y) << endl;
    year=year+1.0/timeSteps;
  }

FUNCTION obs_pred

  ofstream ofs_op("obs_pred_results.dat");
  ofs_op << "survey year sex a_r s_c snum obs pred" << endl;
  year=fyear;

  for(y=startIndex; y<=endIndex; y++) {

    //total observed and predicted catch
    ofs_op << "0 "<< year << " t a c 0 " << TC_obs(y) << " " << TC(y) <<  endl;

    //adult surveys
    for (i=1; i<=numAdSurv; i++)
    ofs_op << i << " "<< year << " 0 a s 0 " << ad_survey_obs(i,y) << " " << ad_survey_est(i,y) << endl;

    //recruit surveys
    for (i=1; i<=numRecSurv; i++)
    ofs_op  << i << " "<< year << " 0 r s 0 " << re_survey_obs(i,y) << " " << re_survey_est(i,y) << endl;

    if (y==startIndex) ofs_op << "0 "<< year << " r r r 0 " <<  R(y) << " " << "NA"  <<  endl;
    else {
      if (SRSwitch==1) ofs_op << "0 "<< year << " r r r 0 " << R(y)<< " "   << SP(y-1)/(SP(y-1)*beta+alpha)  <<  endl;
      if (SRSwitch==2) ofs_op << "0 "<< year << " r r r 0 " << R(y)<< " "   << alpha*SP(y-1)*exp(-beta*(SP(y-1)))  <<  endl;
    } //recruitment deviations

    year=year+1.0/timeSteps;
  }

FUNCTION HPD_estimates
  ofstream ofs_hpd("HPD_results.dat");
  ofs_hpd << "year Adult Spawners Rec RecSurvey1 TC recM adM  F FMSYRatio NMSYRatio FFLimRatio NNLimRatio u0 u1 uAll SREnv MEnvRec MEnvAd" << endl;

  year=fyear; //for outputting the year if multiple time steps per year

  for(y=startIndex;y<=endIndex+projSteps;y++){
    ofs_hpd << year << " " << N(y) << " "<< SP(y)<< " " << R(y) << " " <<R(y)*exp(-sr_time(1)*Mt(1,y)) << " " << TC(y) << " " << Mt(1,y) << " " << Mt(2,y) << " " << F(y)<< " " << F(y)/FMSY<< " " << N(y)/NMSY<< " "<< F(y)/FLim<< " " << N(y)/NLim<< " "<< (sel(1)*F(y))/(sel(1)*F(y)+Mt(1,y))*(1.-exp(-(sel(1)*F(y)+Mt(1,y)))) << " "<< (sel(2)*F(y))/(sel(2)*F(y)+Mt(2,y))*(1.-exp(-(sel(2)*F(y)+Mt(2,y)))) << " "<< u(y) << " " << env(envRecTS,y-envRecLag) << " "<< env(envMTS(1),y-envMLag(1)) << " " << env(envMTS(2),y-envMLag(2)) << endl;
    year=year+1.0/timeSteps;
  }

FUNCTION MSY_estimates
  ofstream ofs_msy("MSY_results.dat");
  {
    //Column headings
    ofs_msy << "Fval\t" << "C_eq\t" << "N_eq\t" << "R_eq\t" << "YPR\t" << "SPR\t" << "SPRatio\t" << "u0_eq\t" << "u1_eq\t" << "uAll_eq\t" << endl;

    Fval=Fval_init;
    for(i=1; i<=Fval_num; i++) {
      ofs_msy << Fval << "\t" << C_eq(i) << "\t" << N_eq(i) << "\t" << R_eq(i) << "\t" <<YPR(i) << "\t" << SPR(i) << "\t" << SPR(i)/SPR0 << "\t" << u0_eq(i) << "\t" << u1_eq(i)<< "\t" << uAll_eq(i) << endl;
      Fval+=Fval_inc;
    }
  }

FUNCTION general_report
  ofstream ofs_gen("gen_results.dat");
  {
    ofs_gen << "Name Value" << endl;
    ofs_gen << "negLL " << negLL <<endl;
    ofs_gen << "Lsa "<< Lsa <<endl;
    ofs_gen << "Lsr "<< Lsr << endl;
    ofs_gen << "Lc "<<  Lc << endl;
    ofs_gen << "Lz "<<  Lz << endl;
    ofs_gen << "Lrdev "<< Lrdev << endl;
    ofs_gen << "Leff "<< Leff << endl;
    ofs_gen << "init_N "<< exp(log_init_N) << endl;
    ofs_gen << "init_R " << exp(log_init_R) << endl;
    for (i=1; i<=numAdSurv; i++) ofs_gen << "qa_i"<<i <<" " << qa(i) << endl;
    for (i=1; i<=numRecSurv; i++) ofs_gen << "qr_i"<<i <<" " << qr(i) << endl;
    ofs_gen << "F_q "<< exp(log_F_q) << endl;
    ofs_gen << "rec_cv "<< exp(log_rec_cv) << endl;
    ofs_gen << "p_rec "<< p_rec << endl;
    ofs_gen << "p_under " << p_under << endl;

    ofs_gen << "SRType " << SRSwitch << endl;
    ofs_gen << "S0 " << exp(log_S0) << endl;
    ofs_gen << "steepness " << exp(log_steep) << endl;
    ofs_gen << "alpha " << alpha << endl;
    ofs_gen << "beta " << beta << endl;
    ofs_gen << "sr_beta_env " << sr_beta_env << endl;
    ofs_gen << "M_beta_env_1 " << M_beta_env(1) << endl;
    ofs_gen << "M_beta_env_2 " << M_beta_env(2) << endl;
    ofs_gen << "sel_1 " << exp(log_sel(1)) << endl;
    ofs_gen << "sel_2 " << exp(log_sel(2)) << endl;
    ofs_gen << "Mr " << Myr(1) << endl;
    ofs_gen << "Ma " << Myr(2) << endl;
    ofs_gen << "Ma " << Myr(2) << endl;
    ofs_gen << "sp_time " << sp_time <<endl;
    ofs_gen  << "MSY " << MSY << endl;
    ofs_gen  << "FMSY " << FMSY << endl;
    ofs_gen  << "FMSYRatio " << FMSYRatio << endl;
    ofs_gen  << "NMSY " << NMSY << endl;
    ofs_gen  << "NMSYRatio " << NMSYRatio << endl;
    ofs_gen  << "RMSY " << RMSY << endl;
    ofs_gen  << "u0MSY " << u0MSY << endl;
    ofs_gen  << "u1MSY " << u1MSY<< endl;
    ofs_gen  << "uMSY " << uMSY<< endl;
    ofs_gen  << "UMSYRatio " << UMSYRatio << endl;
    ofs_gen  << "FLim " << FLim << endl;
    ofs_gen  << "FFLimRatio " << FFLimRatio << endl;
    ofs_gen  << "NLim " << NLim << endl;
    ofs_gen  << "NNLimRatio " << NNLimRatio << endl;
    ofs_gen  << "cLim " << cLim << endl;
    ofs_gen  << "OFL " << OFL<< endl;
    ofs_gen  << "projYears " << projYears<< endl;

    for(ispr=1; ispr<=nspr; ispr++) {
      ofs_gen  << "F"<<SPR_targ(ispr) << "% " << FSPR_ref(ispr)<< endl;
    }
  }

//###################################################################
//###################################################################
REPORT_SECTION
//Call reporting functions
  calculate_tSPR();
  obs_pred();
  MSY_estimates();
  HPD_estimates();
  general_report();

  report << "Likelihood Components" <<endl;
  report << "negLL\t" << negLL <<endl;
  report << "Lsa\t" << Lsa <<endl;
  report << "Lsr\t" << Lsr << endl;
  report << "Lc\t" <<  Lc << endl;
  report << "Leff\t" <<  Leff << endl;
  report << "Lz\t" <<  Lz << endl;
  report << "Lrdev\t" << Lrdev << endl;
  report << "F_pen\t" << F_pen << endl;
  report << "\nParameter Estimates (NOT log space unless marked)" <<endl;
  report << "init_N\t" << exp(log_init_N) << endl;
  report << "init_R\t" << exp(log_init_R) << endl;
  report << "F\t" << F << endl;
  report << "M_rec\t" << Mt(1) << endl;
  report << "M_ad\t" << Mt(2) << endl;
  report << "AveF\t" << mean(F(startIndex,endIndex-1)) << endl;
  report << "AveZ\t" << mean(Z(startIndex,endIndex-1)) << endl;
  report << "AveU\t" << mean(u(startIndex,endIndex)) << endl;
  report << "F_q "<< exp(log_F_q) << endl;
  report << "log_F_dev\t" << log_F_dev << endl;
  report << "mean(log_F_dev)\t" << mean(log_F_dev) << endl;
  report << "rec_dev\t" << exp(log_rec_dev) << endl;
  report << "mean(log_rec_dev)\t" << mean(log_rec_dev) << endl;
  report << "mean(rec_dev)\t" << mean(exp(log_rec_dev)) << endl;
  report << "rec_cv\t" << exp(log_rec_cv) << endl;
  for (i=1; i<=numAdSurv; i++) report << "qa_i"<<i <<"\t" << qa(i) << endl;
  for (i=1; i<=numRecSurv; i++) report << "qr_i"<<i <<"\t" << qr(i) << endl;
  if (SRSwitch==1)	 //Beverton-Holt
  report << "SR=Beverton-Holt\t" << endl;
  if (SRSwitch==2)	 //Ricker
  report << "SR=Ricker\t" << endl;
  report << "S0\t" << exp(log_S0) << endl;
  report << "steepness\t" << exp(log_steep) << endl;
  report << "alpha\t" << alpha << endl;
  report << "beta\t" << beta << endl;
  report << "sr_beta_env\t" << sr_beta_env << endl;
  report << "M_beta_env\t" << M_beta_env << endl;
  report << "M_beta_env_rec\t" << M_beta_env(1) << endl;
  report << "M_beta_env_ad\t" << M_beta_env(2) << endl;
  report << "sel\t" << exp(log_sel) << endl;
  report << "p_rec\t" << p_rec << endl;
  report << "p_under\t" << p_under << endl;
  report << "M\t" << Myr << endl;
  report << "rec_cv\t" << exp(log_rec_cv) << endl;
  report << "sp_time\t" << sp_time <<endl;

  report << "\nReference Point Calculations" <<endl;
  report << "negLL\t" << negLL <<endl;
  report  << "MSY\t" << MSY << endl;
  report  << "uMSY\t" << uMSY<< endl;
  report  << "NMSY " << NMSY << endl;
  report  << "UMSYRatio " << UMSYRatio << endl;
  report  << "NMSYRatio " << NMSYRatio << endl;
  report  << "FMSY\t" << FMSY << endl;
  report  << "FMSYRatio " << FMSYRatio << endl;
  report  << "RMSY\t" << RMSY << endl;
  report  << "u0MSY\t" << u0MSY << endl;
  report  << "u1MSY\t" << u1MSY<< endl;
  report  << "FLim\t" << FLim << endl;
  report  << "FFLimRatio\t" << FFLimRatio << endl;
  report  << "NLim\t" << NLim << endl;
  report  << "NNLimRatio\t" << NNLimRatio << endl;
  report  << "cLim\t" << cLim << endl;
  report  << "OFL\t" << OFL << endl;
  report  << "SPRCurr\t" << SPRCurr<< endl;

  report <<"\n"<< negLL <<endl;
  report  <<MSY << endl;
  report  <<FMSY<< endl;
  report  <<NMSY << endl;
  report  <<FFLimRatio << endl;
  report  <<NNLimRatio << "\n"<< endl;

  for(ispr=1; ispr<=nspr; ispr++) {
    report  << "F"<<SPR_targ(ispr) << "%\t" << FSPR_ref(ispr)<< endl;
  }

//###################################################################
//###################################################################
GLOBALS_SECTION

  #include "admodel.h"
  //define constant variable
  const double MathPI = 3.141592654; //or using MI_PI
  const double EPS = 1e-30;
  const double MathE = 2.71828183;

//###################################################################
//###################################################################
RUNTIME_SECTION
  maximum_function_evaluations 5000,25000,20000,20000,20000,20000
  convergence_criteria 1.0e-8

  //Leave space below this line
Appendix A.2   ADMB Reference Data File (CMSA.dat)

Note: this file selects the stock data to run

#Data file name (uncomment just one filename to run):
#
CMSA_EastStock.dat
#CMSA_WestStock.dat

Appendix A.3   ADMB Data File for East Stock Base Run (CMSA_EastStock.dat)

##########################################################
##########################################################
##########################################################
# Data sources
##########################################################
##########################################################
##########################################################

#Run in testing mode: runs model at initial values and output some values to console (0=off, 1=on)
0

#first year / last year for the model simulation (should be same as catch)
1989 2011

#Retrospective NumYears
0

#Projection NumYears
0

#################################
#Catch Data
#################################

#first / last year in any region of total catch time series
1986 2011

#Total Commercial Catch (in 1 millions of crabs)
#NOAA	
20.67371454	22.9730814	26.24260698	15.01907907	14.03195872	17.11797094	18.36188837	21.14422848	19.00864186	23.81577034	29.32890872	19.81952673	36.56660493	17.03767284	12.2078978	11.97206516	13.60857137	19.71102984	18.2597561	18.28278544	18.58467809	8.362494842	7.267201994	10.19426815	17.06753423	12.62101973

#Commercial catch CV (same for each year since no data on variability):
0.05

#Flag to include effort time series in calcs (adds negLL component for F-deviations)
0

#Trap Effort  (if don't have an effort time series, set all equal to 1 for the total number of years)
#From FL trip ticket program - # traps pulled
5.163355	5.08823	5.875669	4.621626	4.395232	4.483885	5.274843	6.325681	6.053532	6.861174	7.820214	6.298735	8.336047	5.406051	4.409944	4.779467	5.031937	6.07831	5.273512	4.999432	4.992762	3.379471	3.094001	3.260352	4.252351	3.438569

#################################
#Survey Data
#################################

###Adult surveys AND CVs###
#Number of adult surveys 
1

#first / last year in adult surveys
#Note: if catch are different lengths of time, use -999. for missing values
#	Therefore, this is min and max year for any data
1989 2011 

#Standardized (x/mu)

#Glf All Ad, 1990-present (otter trawls 1990-2011, 183m haul seines 1996-2011)
2.801450204	1.09128937	0.793509531	2.132091602	0.640488149	0.994606188	0.597634102	1.452800261	0.460792438	1.916935486	0.788916559	0.413885163	0.277185746	0.390012294	0.671310537	1.141437488	1.0203978	2.297395155	0.738431001	0.186284873	0.235016851	1.173701304	0.7844279
#Glf All Ad CV	
0.374387019	0.40079575	0.406453272	0.339944414	0.445229963	0.50854749	0.334936697	0.144880941	0.166922492	0.116233374	0.121979485	0.15757089	0.171345178	0.152041413	0.127601251	0.10337698	0.106676832	0.094544103	0.114987264	0.207481909	0.177980827	0.110056199	0.118278162

#Adult survey time(s) 
0.0

###recruitment surveys AND CVs###

#Number of recrtui surveys
1

#first /last year in recruit surveys 
#Note: if catch are different lengths of time, use -999. for missing values
#	Therefore, this is min and max year for any data
1989 2011 

#Standardized (x/mu)
#Do all surveys as rows first, then all CVs as rows 2nd
#GlfAll Rec	
1.001005116	0.283812494	1.143517535	0.928502926	0.792187488	0.639873842	0.83191239	0.706005627	0.965528699	2.013047475	1.304789755	1.075983431	0.765711243	0.643204199	1.776157051	1.114431304	0.853290469	1.596605041	0.829762155	0.569516693	0.727422957	1.666886277	0.770845833
#GlfAll Rec CV	
0.366562411	0.440732341	0.288289062	0.297553426	0.277851454	0.242664878	0.130397553	0.115172137	0.11952129	0.109251156	0.113825168	0.107998652	0.112934582	0.113641943	0.091612603	0.093220916	0.09356952	0.091690735	0.101344351	0.106715434	0.1061701	0.098185682	0.149213969

#Recruit survey time
#NOTE: assummed that recruits are surveyed before vulnerable to fishery
#E.g., re_survey_est(r,i)=qr(r,i)*(R(r)*mfexp(-(sr_time(r,i)*Mt(1)));
0.5

#################################
#Fishery params
#################################
#Proportion of recreational harvest per region
.05

#Proportion under reporting per region
#could alternative add this directly to the catch data for month and region specific 
0

#Max F
4

#Max M
4

#Ave Z prior
# calculated as sum of size freq across all years for 183m haul seine and 6.1m otter trawl,
#   using eq 4.4.5.3 in Sparre and venema 1998
# based on fit to simulated data from IBM molt-process temperature
#   for both sexes combined
3.10599543

#+/- 20%
#2.174196801
#4.037794059

#Ave Z CV
.05

#################################
#Life History params
#################################

#Sex ratio (leave at 1)
.5

#Natural mortality per stage 
#3yr mortality 
1.28 .87
#2yr mortality
#1.78 1.22
#3yr constant
#1 1

#Proportion spawning per stage (1st = Recruits, 2nd = Adults)
1
1

#Proportion of the time step before spawning occurs (0=start of year, 1=end of year)
1

#SR formulation (Bev Holt=1, Ricker=2)
2

#################################
#Environmental time series params/data
#################################

#Number of environmental time series
1

#first / last year for the environmental time series
1980 2011

#GlfStream	
-0.982937731	-0.858425202	1.539823242	0.916872095	-0.625035665	0.179729529	0.603164622	-0.126476777	-0.255735071	0.374830432	-0.279651943	0.094259549	0.714082581	-0.603017897	1.486887988	1.209786676	-0.394459622	2.378390258	-0.262426043	-1.043837322	-1.47568612	-0.884639379	0.532935588	0.900764451	1.256223547	0.468701733	-1.223624997	-1.287573367	-0.551710175	0.794235242	-1.015937809	-1.579512414
#AtlStream	-1.371195991	-0.934229209	2.01234639	1.041339067	-0.036041279	0.306438513	-0.024093324	-0.441231402	-0.770842701	-1.189193943	-0.538279926	0.088980559	0.67933287	-0.709307014	1.752663352	1.075693309	-0.201937715	2.291092663	-0.5932093	0.204586106	-1.388224535	0.105886731	0.659193015	-0.522463434	0.849843471	1.2054167	-1.054161503	-1.176871965	0.228422422	0.456717828	-1.158008571	-0.848661184
#FLPrecip	-1.676902843	0.69198606	1.344133542	0.252757366	-1.484070733	0.471032601	0.544683754	-0.693994728	-0.229322908	-0.75961121	0.908922184	-0.003013002	-0.072646819	-0.040508134	1.303960186	1.234326369	-0.148976196	0.919635079	0.116167955	-1.33141198	-0.904235292	0.349173421	1.590530127	-0.455632815	2.051184611	-1.035468256	-1.658155277	-0.336451858	0.401398784	0.310339177	-1.78938824	0.129559074

#Environment series CV 
.1

#Time series which influences recruitment:
#Use 3, assumming the num of eggs are survival of early recruits are due to rain 6 months around spawn peak
1

#Lag in environment influence on recruitment (# time steps)
1

#Time series that influences mortality (one for each stage):
1 1

#Lag in enviro influence on mortality (one for each stage):
0 0

#################################
#Projections time series
#################################

#Environmental time series annomolies for projection years +1 (+1 is for recruit calc) 
#Note: must be same number of series as in environment section above with mean=0 (average)
#0 0 0 0 0 0 0 0 0 0 0
0
#Effort Deviation from year before terminal yr (0) to FMSY (1)
#E.g., 0 .25 .5 .75 1 would be step increase from year before terminal yr F -> FMSY
#.25 .5 .75 1 1 1 1 1 1 1 1
0

##########################################################
##########################################################
##########################################################
#Parameters and flags
# Format: 
#	1st row: initial parameter estimates vector by sub-stock (or stage, e.g., selectivity)
#  2nd row: min bound, max bound, phase of estimatation
# note: if phase <0, then initial value will be held constant
##########################################################
##########################################################
##########################################################

#################################
#Initial values, bounds, and phase
#################################

#init_N (first row is vector of initial parameter guesses by sub-stock)
5
.01 100 1

#init_R 
#Note: only used if the SR lag (in years) is >0
50
.01 500 1

################
# F params
################
## F=(q*Effort)*exp(Fdev)  where 1st Fdev=0 so q scales to the initial year F
#F q 
1
.00001 10 2

#F_dev
1 
.1 5 5

#effort_cv
.2
.01 1 -1

################
# Recruit params
################

#rec_dev (expected: log(rec_dev)=0 / rec_dev=1)
1 
.1  5 4

#rec_cv 
.5
.3 1 -1

#Stock Recruitment S0
5 
.1 10000 3

#S-R steepness
.8
.2 .9999 3

########
# Environment params
########
#sr_beta_env (environmental link parameter for recruitment: R=R*exp(sr_beta_env*env) )
0
-20 20 -3

#M_beta_env (environmental link parameter for yearly M: Mt=M*exp(M_beta_env*env) )
#note: one for each stage
0 0
-20 20 -3

#M_cv 
.1
.3 1 -1

#sel (vulnerability)
#NOTE: currently the modle only accomodates surveying the recruits prior to being vulnerable
# if this is not the case (i.e., sel > (1-sr_time)), need to re-code the index calcs
0.3 1.0
0.1 1 -1

#################################
#Likelihood weights
#################################

#Landings weight(s) lambda
1.0
#Adult survey weight(s) lambda
1.0
#Recruit survey weight(s) lambda 
1.0
#recruitment deviation weight(s) lambda 
1.0
#effort residuals weight(s) lambda 
1.0
#Z prior weight(s) lambda 
1.0

#################################
#Additional param control flags not addressed in data section
#################################

#Bias correction adjustment for predicted recruitment: biasAdj*(0.5*var)
#can turn off by setting=0 or turn on to whatever proportion by setting =1
1

#################################
#Reference point calcs 
#################################

#variables to control F for females used in reference point calculations
#FSPR_init FSPR_max FSPR_inc
0 6.0 0.01

#SPR targets for calculating F reference points
#number of SPR targets
5

#targets
0.05 0.1 0.2 0.3 0.4 

#################################
#EOF I/O test
#################################
#EOF number
12345

Appendix A.4   ADMB Data File for West Stock Base Run (CMSA_WestStock.dat)

##########################################################
##########################################################
##########################################################
# Data sources
##########################################################
##########################################################
##########################################################

#Run in testing mode: runs model at initial values and output some values to console (0=off, 1=on)
0
 
#first year / last year for the model simulation (should be same as catch)
1985 2011

#Retrospective NumYears
0

#Projection NumYears
0

#################################
#Catch Data
#################################

#first / last year in any region of total catch time series
1985 2011

#Total Commercial Catch (in 1 millions of crabs)
#rows=spatial regions, columns=time steps (steps/year * years)
#NOAA	
#Grand Total	
98.45087442	121.6449721	159.966614	142.6423512	114.9461349	133.5054744	135.1838977	144.407093	111.3985395	111.178993	107.7819372	117.0608488	126.4343372	140.4308372	157.5395901	118.6543769	129.1933266	143.8513492	127.9418075	113.460141	115.9262415	126.9902281	114.5286733	127.0669509	108.7038416	89.52476457	89.35450346
#AL	5.435304654	6.183193021	7.115025582	9.647239535	8.342374421	8.282074419	5.924637207	7.307462788	6.685593024	5.477879064	6.590725582	7.937344183	7.04770465	9.082044183	11.17805116	8.583411628	5.083004651	6.036727904	7.704506973	5.383213957	3.709179072	6.479458138	5.069918609	4.213295355	2.554772091	3.030869774	3.545802329
#LA	69.08582093	91.41135349	122.8216488	105.7652395	86.10353721	106.1285116	114.8085	118.7742884	88.97858139	92.92577906	84.88989534	93.86322558	102.0428465	111.1367279	131.4285739	100.0005699	108.5660359	120.5068515	108.7324215	97.57581307	103.8837973	114.6192514	98.28634307	115.5828393	99.48527651	77.48651108	82.15458485
#MS	1.110269768	3.367879068	2.423293026	1.919886048	1.210693022	1.309462792	1.017881395	0.707223256	0.482860465	0.444767442	1.031639535	0.940999999	1.623460464	1.884706976	2.006048838	1.424160466	1.196513953	2.008011627	2.216918607	1.757376743	1.558572093	1.903290698	1.469530232	1.275102327	0.943355813	0.882525581	1.079120931
#TX/WestSubregion	22.81947907	20.68254651	27.60664652	25.30998604	19.28953023	17.78542558	13.43287907	17.6181186	15.25150465	12.33056744	15.26967675	14.31927907	15.72032558	18.32735814	12.92691628	8.646234885	14.3477721	15.29975814	9.287960463	8.743737209	6.77469302	3.988227902	9.702881396	5.995713951	5.720437205	8.124858144	2.574995349
#CentralSubregion	75.63139535	100.9624256	132.3599674	117.3323651	95.65660465	115.7200488	121.7510186	126.7889744	96.14703488	98.84842557	92.51226046	102.7415698	110.7140116	122.1034791	144.6126739	110.008142	114.8455545	128.551591	118.6538471	104.7164038	109.1515485	123.0020002	104.8257919	121.071237	102.9834044	81.39990644	86.77950811

#Commercial catch CV (same for each year since no data on variability):
0.05

#Flag to include effort time series in calcs (adds negLL component for F-deviations)
0

#Trap Effort  (if don't have an effort time series, set all equal to 1 for the total number of years)
#
1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1

#################################
#Survey Data
#################################

###Adult surveys AND CVs###
#Number of adult surveys 
1

#first / last year in adult surveys
#Note: if catch are different lengths of time, use -999. for missing values
#	Therefore, this is min and max year for any data
1985 2011 

#Standardized (x/mu)
#Do all surveys as rows first, then all CVs as rows 2nd
#All Adults	mean	
1.766147166	1.768107473	1.596853607	1.10483864	0.999790093	1.526579353	1.42450192	1.039353467	1.466427913	0.939851755	0.436654142	0.613031703	0.786854207	0.964869567	0.690056536	0.641226945	0.634253061	0.78304312	0.753201752	0.786065262	0.740924042	1.358614501	0.97159447	0.728005718	0.890133492	0.786544662	0.802475432
0.064247085	0.065432786	0.061646786	0.064997891	0.068393927	0.065248206	0.065254712	0.06931112	0.065986747	0.070719855	0.082478639	0.077020954	0.070775011	0.068675335	0.073868918	0.076923252	0.073756897	0.070873845	0.071980759	0.06876141	0.076525005	0.064977194	0.068481944	0.072156128	0.068916816	0.075843186	0.073334971
#AL adults	mean	
#1.934161061	1.905566414	0.651811671	0.847748157	1.91460443	1.269750201	1.119446395	1.3313397	1.117875209	0.707451824	0.806645291	0.89329348	0.805844532	0.894965339	0.659924503	0.806577107	0.986254122	0.809719193	1.578323541	1.725442417	0.407125917	0.740273822	0.328870487	0.315010261	0.432551229	1.109012323	0.900411374
#0.240723558	0.209773392	0.255863755	0.263971498	0.253668936	0.24920737	0.24784996	0.241078466	0.255918017	0.287586546	0.273385335	0.309036602	0.278979121	0.274338189	0.350830256	0.32896739	0.209416765	0.244914758	0.214925863	0.191024744	0.307338662	0.268476953	0.368926996	0.381085994	0.331218081	0.270736791	0.286838964
#LS Adults	mean	
#1.965712926	1.844279267	1.34838459	0.894541797	0.901581753	1.714053028	1.626328809	0.741896556	0.773819759	0.581447302	0.271069218	0.409413256	0.605203225	0.656865133	0.710436113	0.756996339	0.468357536	0.802997954	0.490993146	0.597184321	1.209127263	2.559258298	1.154544038	0.735864931	1.191863172	0.848387437	1.139392831
#0.100009391	0.110027442	0.09829643	0.106866881	0.117501811	0.097179463	0.096538353	0.111266888	0.110102397	0.120857132	0.1396424	0.131976617	0.110478769	0.112045926	0.110842472	0.106775877	0.118157724	0.108309045	0.120157497	0.114157905	0.111336732	0.100705248	0.103905552	0.107584485	0.097738162	0.112830586	0.104688292
#MS adults	mean	
#1.062864439	1.479250935	0.79254653	1.23232732	2.504375184	1.345340552	0.22209675	0.526642023	0.689690062	0.852902055	1.009943561	2.257556179	1.944202645	0.82824453	0.762871277	1.01176093	1.605173552	0.707178272	0.634301661	0.670106532	0.416955051	0.816896111	0.84321792	0.262125511	0.763117353	1.480115809	0.278197256
#0.548760991	0.410409083	0.502978494	0.4184149	0.377290973	0.706681493	1.74194089	2.30971521	0.825390596	0.684400629	0.873328024	0.535363604	0.576940588	0.902191271	0.929316607	0.708361791	0.60228561	0.924976488	0.853158878	0.827994075	1.51691707	1.62494303	0.779682238	2.0457767	2.2257102	0.642529604	2.07795042
#TX Adults	mean	
#1.60360775	1.789673715	2.030992469	1.305510919	0.926925935	1.388987099	1.298031926	1.285060252	2.184216893	1.266426715	0.531273635	0.717301929	0.904683137	1.24191986	0.625789687	0.484048002	0.69883636	0.738857404	0.846177418	0.770935246	0.424828449	0.424765162	0.923334521	0.764766034	0.655271569	0.670664734	0.497113178
#0.09185961	0.083614965	0.083570045	0.086875809	0.091534378	0.08856159	0.087259405	0.093678848	0.0844765	0.090921221	0.113393377	0.10103429	0.094411799	0.088617666	0.103466211	0.112489994	0.103806482	0.098546766	0.095217299	0.097702481	0.116794742	0.117793316	0.094947395	0.098760169	0.104670116	0.104972032	0.112319954
#Central Adults	mean	
#1.933002267	1.812643788	1.207798908	0.90395831	1.083142796	1.666135116	1.561080499	0.817750948	0.824945217	0.61331737	0.347652511	0.49286234	0.660189478	0.688459766	0.738852761	0.772606991	0.563642355	0.798951724	0.638545961	0.777571467	1.067099317	2.264010375	1.016007994	0.678738642	1.082610508	0.88827046	1.100152133
#0.092088354	0.096776293	0.090685224	0.096700229	0.101655466	0.091964072	0.092193626	0.102222007	0.10255344	0.110612017	0.121496563	0.117061245	0.102588433	0.103474237	0.106829555	0.101514333	0.10298717	0.098932384	0.105915274	0.098676099	0.104150861	0.089676908	0.09834757	0.103682887	0.095195677	0.103757257	0.09826679

#Adult survey time(s) 
0.0

###recruitment surveys AND CVs###

#Number of recruit surveys
1

#first /last year in recruit surveys 
#Note: if catch are different lengths of time, use -999. for missing values
#	Therefore, this is min and max year for any data
1985 2011 

#Standardized (x/mu)
#Do all surveys as rows first, then all CVs as rows 2nd
#All Recruits	mean	
1.281060449	0.961905987	1.684921265	1.366795752	1.546392544	1.888464386	1.267444476	1.149981792	1.517864986	1.233561953	0.993402856	0.881826873	1.393501748	0.659525956	1.128206135	0.641399509	0.55726424	0.489175454	0.722851228	0.794348598	0.893163237	0.599717045	0.613062775	0.716521884	0.371561463	0.671836817	0.974240592
0.082604791	0.085677496	0.082567205	0.085648035	0.087637115	0.086046495	0.085287309	0.087673689	0.085091555	0.085467383	0.089105278	0.087590841	0.080807254	0.083848689	0.087980306	0.092260579	0.088100563	0.087623753	0.084544347	0.082349279	0.085764847	0.088823729	0.087530303	0.086073428	0.094116767	0.095835502	0.118244586
#AL Recruits	mean	
#2.259304668	1.362175172	2.293538512	1.228282876	1.047260519	2.038404855	0.875922401	0.508953165	1.117391961	1.162530918	0.526935763	1.281291823	1.694670196	0.643649025	0.836917276	0.272098583	0.507071476	0.268250578	1.897288978	0.435408108	0.686091597	0.930567871	0.539312174	0.24214521	0.236441482	0.743370888	1.364723925
#0.27103365	0.257760793	0.22508085	0.341360378	0.309843555	0.278263529	0.282205007	0.300645097	0.315732484	0.293787719	0.328172826	0.273370631	0.275979873	0.380419755	0.373862642	0.349074225	0.324320952	0.348294607	0.255216171	0.283016542	0.283983707	0.285358928	0.310096529	0.369042991	0.331600568	0.346818871	0.453133848
#LS Recruits	mean	
#1.021414525	1.089342144	1.446027319	0.900833691	1.169274102	1.771677384	1.046041028	0.77154291	1.617251266	1.172011336	1.063405894	0.880954991	1.684168453	0.492378093	1.30873406	0.878713542	0.651849656	0.48353943	0.766854474	1.15749317	1.14832627	0.718343819	0.846465708	0.738274786	0.48622214	0.648777584	1.040082225
#0.148007688	0.166945921	0.166949606	0.145134448	0.154478846	0.147839232	0.146847808	0.15116891	0.151881965	0.1429286	0.14093523	0.152059651	0.136914614	0.126856454	0.132646646	0.133938077	0.12807167	0.129960964	0.134588077	0.12312838	0.132417641	0.128362326	0.12625073	0.128244464	0.125488848	0.150874589	0.145392029
#MS Recruits	mean	
#1.744040989	1.488430434	1.384528336	1.303839497	2.199279331	2.999370682	1.06599792	0.437624839	0.4347188	1.701416109	1.736085156	1.646329799	1.799698261	0.423138335	0.808186704	0.277205678	0.722019716	0.718267688	0.66186235	0.27862969	1.423439987	0.320311528	0.732136403	0.238224258	0.079569454	0.116903794	0.258744261
#0.421007152	0.411593799	0.414375592	0.409850175	0.669082108	0.664858803	0.729055698	0.941717094	0.917099323	0.611029607	0.712901043	0.600551221	0.676986612	0.85396163	1.24038984	1.53801629	0.827982523	0.679869851	0.629009647	5.10481154	0.654018038	1.4796399	0.773689055	1.57703416	2.26935799	1.43372723	1.61123301
#TX Recruits	mean	
#1.329779181	0.799638554	1.718275586	1.891654321	1.792226945	1.894603557	1.715090415	1.717031845	1.466155363	1.297676432	0.973281401	0.743362651	1.029551594	1.001125473	0.996484034	0.502342569	0.519036216	0.60052039	0.551143019	0.521986637	0.577776193	0.416246246	0.390774299	0.827163645	0.257903773	0.76197142	0.707198241
#0.099693156	0.10958801	0.095740566	0.102084033	0.096841744	0.104094755	0.098776698	0.101428589	0.10031809	0.102181046	0.107205872	0.114912105	0.105115803	0.100718735	0.114655565	0.127421231	0.123059898	0.111299474	0.122000902	0.113793704	0.12364692	0.132919422	0.127831018	0.116281629	0.155926345	0.120774904	0.247378896
#Central Recruits	mean	
#1.206682086	1.144075042	1.59823793	0.951455877	1.191727807	1.845703122	1.015123786	0.708258572	1.470794626	1.178639796	1.004934335	0.954383164	1.690392153	0.496441075	1.30751928	0.789672838	0.627992222	0.446685439	0.87599117	1.013097585	1.122172862	0.73213531	0.813587297	0.66405557	0.424830258	0.638336024	1.087074777
#0.126620075	0.13393635	0.127105647	0.129497784	0.133404695	0.127576894	0.128385577	0.133644844	0.12971474	0.126359249	0.129677474	0.13004819	0.119709539	0.120068466	0.120271754	0.12041334	0.118559533	0.122629509	0.118822676	0.110880573	0.116603728	0.114173279	0.115956484	0.115138876	0.117095314	0.138605109	0.139397655

#Recruit survey time
#NOTE: assummed that recruits are surveyed before vulnerable to fishery
#E.g., re_survey_est(r,i)=qr(r,i)*(R(r)*mfexp(-(sr_time(r,i)*Mt(1)));
0.5

#################################
#Fishery params
#################################
#Proportion of recreational harvest per region
.05

#Proportion under reporting per region
#could alternative add this directly to the catch data for month and region specific 
0

#Max F
4

#Max M
4

#Ave Z prior
# calculated as sum of size freq across all years for 183m haul seine and 6.1m otter trawl,
#   using eq 4.4.5.3 in Sparre and venema 1998
# based on fit to simulated data from IBM molt-process temperature
#   for both sexes combined
2.06290806

#+/- 20%
#1.444035642
#2.681780478

#Ave Z CV
.05

#################################
#Life History params
#################################

#Sex ratio (leave at 1)
.5

#Natural mortality per stage 
#3yr mortality
1.46 0.77
#2yr mortality
#1.95 1.05
#3yr constant mortality
#1 1

#Proportion spawning per stage (1st = Recruits, 2nd = Adults)
1
1

#Proportion of the time step before spawning occurs (0=start of year, 1=end of year)
1

#SR formulation (Bev Holt=1, Ricker=2)
2

#################################
#Environmental time series params/data
#################################

#Number of environmental time series
1

#first / last year for the environmental time series
1980 2011

#These are for adjusted year, not calendar year (July->June)
#TXPrecip	0.222751937	0.426310576	-0.652901173	-1.144249611	0.519315816	0.143785224	1.500257877	-1.079321424	-0.001864492	-0.60727596	0.670229979	2.411007303	0.170107462	-0.803815336	0.503522473	-1.553121704	1.681003909	-0.779247914	0.772009298	-1.005619159	-0.198403867	-0.501987009	0.594772898	0.891336776	0.326286072	-1.270596352	1.5581668	-0.428284744	-0.700281201	0.84746638	-1.556631336	-0.954729499
#LAPrecip	-1.649545417	-1.182987328	2.184253182	-0.287871355	-0.334316052	0.009796928	0.463688282	-0.407149781	1.187803326	0.311687457	1.681806009	0.542855379	0.894357288	-0.585539638	0.598800127	-0.916985883	1.262748178	-0.189704155	0.473188334	-1.767768281	0.148075457	-0.048258943	1.006246784	1.130803016	-0.534872697	-1.353988256	0.179742296	-0.4229832	-0.640428825	0.54391094	-2.040103094	-0.257260078
#ALPrecip	-1.391722302	-0.107246552	1.630768315	-0.027794444	-0.874180097	-1.046326332	0.356224079	-1.878366465	1.211437744	1.327305402	0.586855893	-0.792520986	0.526163311	0.105729238	0.147662295	0.813073701	0.939976374	0.525059809	0.287806986	-1.71284124	-0.026690943	-0.711965376	1.64180333	-0.038829459	0.94439038	-0.588373208	-1.530763492	-0.655686799	0.640927467	1.48731312	-1.105915413	-0.683274337
#MSPrecip	-1.432329829	-0.864596177	2.769103051	-0.181684961	-0.635260321	-0.728013934	0.370759634	-1.599490186	1.261602025	0.973148482	1.520496725	-0.499697348	0.321834651	0.186271679	0.273928939	-0.326421368	1.010861489	-0.182704232	0.043573813	-1.56381572	-0.322344286	0.464532517	1.559229002	0.639847038	-0.126644356	-0.730052475	-1.173435129	-0.447714554	0.398279936	0.748908978	-1.337537675	-0.390635408
#WestPrecip	-1.395710381	-0.87110158	1.80866847	-0.506862226	-0.230291985	-0.07259123	0.759587832	-0.793846979	1.071281146	0.243041627	1.588571371	0.90657151	0.780987342	-0.621123762	0.606726586	-1.044038179	1.495057617	-0.310402757	0.567398784	-1.807325494	0.043010926	-0.169393174	1.120751788	1.120089603	-0.282425741	-1.426624491	0.333825233	-0.497343014	-0.590969438	0.761895721	-2.094163042	-0.493252083
#centralSubPrecip(LA,AL,MS)	-1.677173766	-1.125848061	2.273670734	-0.271847519	-0.409562225	-0.123619227	0.466442223	-0.621971142	1.240770462	0.448534606	1.654655135	0.385919284	0.857351226	-0.497838981	0.56382881	-0.781208528	1.268167936	-0.144869475	0.444394671	-1.815595503	0.104407341	-0.057936253	1.133815232	1.051417119	-0.416786265	-1.30190849	-0.042416946	-0.45790163	-0.491426897	0.648800001	-1.996006399	-0.308257474
#westStream	
-1.347464185	-0.082128462	0.370539104	-0.267880573	-0.59635016	-0.675992103	0.528157922	-0.501676926	-1.023916844	0.757663959	0.320985794	2.53148789	0.747354456	0.022446003	0.822725918	-0.912776692	0.421701302	1.106280384	0.65159242	-1.804756843	-0.774402102	0.694138381	1.310355869	0.878248549	1.355458641	-1.387474055	-0.691620993	0.259773332	-0.827719572	0.674032765	-1.43384874	-1.124934438
#westSubStream(TX)	-1.026546614	0.547815074	-0.406801755	-0.505120733	-0.555864394	-0.834193296	0.262860576	-0.30828085	-0.839352363	0.049416073	0.329000755	3.391508469	0.648105836	-0.251591495	0.87952503	-1.094468658	0.198289033	0.248161334	0.918544526	-1.214136636	-0.145692527	0.765279979	1.030129782	0.410535577	1.199249487	-1.471811055	0.058375104	1.122777218	-1.13004559	0.119587471	-1.182695121	-1.21256024
#centralSubStream(LA,AL,MS)  -1.285957009	-1.649637591	2.030513586	0.463535429	-0.318126139	0.068515657	0.697193704	-0.623317207	-0.849597059	1.526897843	0.145243638	0.015577809	0.533204293	0.405464709	0.330015972	-0.238104586	0.519194734	1.655179695	0.049346883	-1.794872443	-1.246813314	0.300862681	1.079385525	1.084774939	0.973407425	-0.740561339	-1.385348872	-0.97985674	-0.116039912	1.115943775	-1.202570581	-0.563455505

#Environment series CV 
.1

#Time series which influences recruitment:
1

#Lag in environment influence on recruitment (# time steps)
1

#Time series that influences mortality:
1 1

#Lag in enviro influence on mortality
0 0

#################################
#Projections time series
#################################

#Environmental time series annomolies for projection years +1 (+1 is for recruit calc) 
#Note: must be same number of series as in environment section above with mean=0 (average)
#0 0 0 0 0 0 0 0 0 0 0
0
#Effort Deviation from year before terminal yr (0) to FMSY (1)
#E.g., 0 .25 .5 .75 1 would be step increase from year before terminal yr F -> FMSY
#.25 .5 .75 1 1 1 1 1 1 1 1
0

##########################################################
##########################################################
##########################################################
#Parameters and flags
# Format: 
#	1st row: initial parameter estimates vector by sub-stock (or stage, e.g., selectivity)
#  2nd row: min bound, max bound, phase of estimatation
# note: if phase <0, then initial value will be held constant
##########################################################
##########################################################
##########################################################

#################################
#Initial values, bounds, and phase
#################################

#init_N (first row is vector of initial parameter guesses by sub-stock)
100
.01 500 1

#init_R 
#Note: only used if the SR lag (in years) is >0
50
.01 1000 1

################
# F params
################
## F=(q*Effort)*exp(Fdev)  where 1st Fdev=0 so q scales to the initial year F
#F q 
1
.00001 10 2

#F_dev
1 
.1 5 5

#effort_cv
.2
.01 1 -1

################
# Recruit params
################

#rec_dev (expected: log(rec_dev)=0 / rec_dev=1)
1 
.1  5 4

#rec_cv 
.5
.3 1 -1

#Stock Recruitment S0
100
.1 1000 3

#S-R steepness
.8 
.2 .9999 3

########
# Environment params
########
#sr_beta_env (environmental link parameter for recruitment: R=R*exp(sr_beta_env*env) )
0
-20 20 -3

#M_beta_env (environmental link parameter for yearly M: Mt=M*exp(M_beta_env*env) )
#note: one for each stage
0 0
-20 20 -3

#M_cv 
.1
.3 1 -1

#sel (vulnerability)
#NOTE: currently the modle only accomodates surveying the recruits prior to being vulnerable
# if this is not the case (i.e., sel > (1-sr_time)), need to re-code the index calcs
0.3 1.0 
0.1 1 -1

#################################
#Likelihood weights
#################################

#Landings weight(s) lambda
1.0
#Adult survey weight(s) lambda
1.0
#Recruit survey weight(s) lambda 
1.0
#recruitment deviation weight(s) lambda 
1.0
#effort residuals weight(s) lambda 
1.0
#Z prior weight(s) lambda 
1.0

#################################
#Additional param control flags not addressed in data section
#################################

#Bias correction adjustment for predicted recruitment: biasAdj*(0.5*var)
#can turn off by setting=0 or turn on to whatever proportion by setting =1
1

#################################
#Reference point calcs 
#################################

#variables to control F for females used in reference point calculations
#FSPR_init FSPR_max FSPR_inc
0 6.0 0.01

#SPR targets for calculating F reference points
#number of SPR targets
5

#targets
0.05 0.1 0.2 0.3 0.4 

#################################
#EOF I/O test
#################################
#EOF number
12345

Appendix B.1 Input data and parameters for ASPIC surplus production model base configuration for Eastern GOM Stock

FIT  ## Run type (FIT, BOT, or IRF)
"sp-01"
LOGISTIC  YLD    SSE  
103  ## Verbosity
1000  90  ## Number of bootstrap trials, <= 1000
1  2000000  ## 0=no MC search, 1=search, 2=repeated srch; N trials
1.0000E-08  ## Convergence crit. for simplex
3.0000E-08  20  ## Convergence crit. for restarts, N restarts
1.0000E-04  24  ## Conv. crit. for F; N steps/yr for gen. model
6.0000  ## Maximum F when cond. on yield
1.0  ## Stat weight for B1>K as residual (usually 0 or 1)
1  ## Number of fisheries (data series)
1.0000E+00    ## Statistical weights for data series
0.75  ## B1/K (starting guess, usually 0 to 1)
9.89205E+06  ## MSY (starting guess)
9.89205E+07  ## K (carrying capacity) (starting guess)
4.6081E-08    ## q (starting guesses -- 1 per data series)
1  1  1  1    ## Estimate flags (0 or 1) (B1/K,MSY,K,q1...qn)
9.89205E+05  1.97841E+08  ## Min and max constraints -- MSY
9.89205E+06  1.97841E+09  ## Min and max constraints -- K
3921295  ## Random number seed
62  ## Number of years of data in each series
"Florida Index"
CC
   1950    -2.439024E+00     7.186200E+05
   1951    -2.439024E+00     2.184000E+06
   1952    -2.439024E+00     2.098740E+06
   1953    -2.439024E+00     3.314115E+06
   1954    -2.439024E+00     3.048465E+06
   1955    -2.439024E+00     5.202645E+06
   1956    -2.439024E+00     3.915975E+06
   1957    -2.439024E+00     5.577180E+06
   1958    -2.439024E+00     9.128910E+06
   1959    -2.439024E+00     1.459353E+07
   1960    -2.439024E+00     1.958513E+07
   1961    -2.439024E+00     1.799133E+07
   1962    -2.439024E+00     1.087433E+07
   1963    -2.439024E+00     1.381002E+07
   1964    -2.439024E+00     1.478558E+07
   1965    -2.439024E+00     2.163966E+07
   1966    -2.439024E+00     1.737540E+07
   1967    -2.439024E+00     1.468173E+07
   1968    -2.439024E+00     9.458505E+06
   1969    -2.439024E+00     1.216341E+07
   1970    -2.439024E+00     1.552593E+07
   1971    -2.439024E+00     1.289264E+07
   1972    -2.439024E+00     1.120697E+07
   1973    -2.439024E+00     1.007843E+07
   1974    -2.439024E+00     1.064049E+07
   1975    -2.439024E+00     1.344851E+07
   1976    -2.439024E+00     1.265093E+07
   1977    -2.439024E+00     1.662381E+07
   1978    -2.439024E+00     1.228596E+07
   1979    -2.439024E+00     1.176797E+07
   1980    -2.439024E+00     1.185724E+07
   1981    -2.439024E+00     1.555080E+07
   1982    -2.439024E+00     9.370517E+06
   1983    -2.439024E+00     9.841806E+06
   1984    -2.439024E+00     1.358693E+07
   1985    -2.439024E+00     1.290458E+07
   1986    -2.439024E+00     8.036358E+06
   1987    -2.439024E+00     1.094588E+07
   1988    -2.439024E+00     1.092292E+07
   1989     6.832805E+00     8.607252E+06
   1990     2.661681E+00     7.260622E+06
   1991     1.935389E+00     5.496715E+06
   1992     5.200223E+00     8.036314E+06
   1993     1.562166E+00     8.949551E+06
   1994     2.425869E+00     8.887131E+06
   1995     1.457644E+00     9.219875E+06
   1996     3.543415E+00     1.309866E+07
   1997     1.123884E+00     9.787235E+06
   1998     4.675452E+00     1.350592E+07
   1999     1.924187E+00     1.172794E+07
   2000     1.009476E+00     6.901276E+06
   2001     6.760628E-01     4.878993E+06
   2002     9.512495E-01     5.845369E+06
   2003     1.637343E+00     7.586642E+06
   2004     2.783994E+00     8.487322E+06
   2005     2.488775E+00     7.738503E+06
   2006     5.603403E+00     9.040658E+06
   2007     1.801051E+00     6.415316E+06
   2008     4.543533E-01     2.796251E+06
   2009     5.732118E-01     3.532425E+06
   2010     2.862686E+00     6.046601E+06
   2011     1.913239E+00     7.174866E+06

Appendix B.2 Input data and parameters for ASPIC surplus production model base configuration for Western GOM Stock

FIT  ## Run type (FIT, BOT, or IRF)
"sp-01"
LOGISTIC  YLD    WTDSSE  
103  ## Verbosity
1000  90  ## Number of bootstrap trials, <= 1000
0  2000000  ## 0=no MC search, 1=search, 2=repeated srch; N trials
1.0000E-08  ## Convergence crit. for simplex
3.0000E-08  8  ## Convergence crit. for restarts, N restarts
1.0000E-04  24  ## Conv. crit. for F; N steps/yr for gen. model
6.0000  ## Maximum F when cond. on yield
1.0  ## Stat weight for B1>K as residual (usually 0 or 1)
1  ## Number of fisheries (data series)
1.0000E+00    ## Statistical weights for data series
0.75  ## B1/K (starting guess, usually 0 to 1)
4.64508E+07  ## MSY (starting guess)
4.64508E+08  ## K (carrying capacity) (starting guess)
2.7894E-09    ## q (starting guesses -- 1 per data series)
1  1  1  1    ## Estimate flags (0 or 1) (B1/K,MSY,K,q1...qn)
4.64508E+06  9.29016E+08  ## Min and max constraints -- MSY
4.64508E+07  9.29016E+09  ## Min and max constraints -- K
3921295  ## Random number seed
62  ## Number of years of data in each series
"Western Combined INDEX"
CC
1950	-2.43902439	19420380	0
1951	-2.43902439	14866530	0
1952	-2.43902439	13140750	0
1953	-2.43902439	15443085	0
1954	-2.43902439	13703445	0
1955	-2.43902439	21088725	0
1956	-2.43902439	17469480	0
1957	-2.43902439	19427835	0
1958	-2.43902439	23627415	0
1959	-2.43902439	30842595	0
1960	-2.43902439	37169790	0
1961	-2.43902439	37686390	0
1962	-2.43902439	27557670	0
1963	-2.43902439	28199535	0
1964	-2.43902439	26782665	0
1965	-2.43902439	39086145	0
1966	-2.43902439	32633790	0
1967	-2.43902439	29066205	0
1968	-2.43902439	27345990	0
1969	-2.43902439	35055195	0
1970	-2.43902439	35795025	0
1971	-2.43902439	35340270	0
1972	-2.43902439	37062375	0
1973	-2.43902439	45771810	0
1974	-2.43902439	42471975	0
1975	-2.43902439	40771605	0
1976	-2.43902439	38482185	0
1977	-2.43902439	46780230	0
1978	-2.43902439	40245393.3	0
1979	-2.43902439	45835585.95	0
1980	-2.43902439	45005747.85	0
1981	-2.43902439	44549579.55	0
1982	-2.43902439	38793170.85	0
1983	-2.43902439	42494597.25	0
1984	-2.43902439	59143614.6	0
1985	4.307676015	58644988.5	242.2663769
1986	4.312457251	55665646.05	233.5657594
1987	3.894764895	82389458.55	263.135305
1988	2.69472839	83244848.4	236.7017504
1989	2.438512422	58425216.15	213.7789398
1990	3.723364276	61263699	234.8890922
1991	3.474394927	69124282.5	234.8422569
1992	2.535008456	73283362.95	208.1585081
1993	3.576653446	68840932.65	229.6606352
1994	2.292321354	55904844.45	199.9480978
1995	1.065010102	57094206.75	146.9998962
1996	1.495199276	65535013.95	168.5707476
1997	1.919156602	67307339.4	199.6365743
1998	2.353340407	70958647.2	212.0305436
1999	1.683064722	72445591.05	183.2637313
2000	1.563968159	72342378.15	168.9992307
2001	1.546958685	57225395.85	183.8208311
2002	1.909861268	69319649.7	199.0801737
2003	1.837077444	67158598.5	193.0043763
2004	1.917232346	63609825.3	211.5000395
2005	1.80713181	52543234.8	170.7627998
2006	3.313693905	70855203.3	236.8525663
2007	2.36974261	60862579.05	213.2297706
2008	1.775623702	51722924.4	192.0673567
2009	2.171057298	64335784.8	210.5472585
2010	1.918401615	43301931.75	173.8468649
2011	1.957257151	58385875.8	185.9421083